Journal  of  Power  Sources  198  (2012)  329-337 


ELSEVIER 


Contents  lists  available  at  SciVerse  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


in 

SllbaidtS 


Simplification  and  order  reduction  of  lithium-ion  battery  model  based  on 
porous-electrode  theory 

Thanh-Son  Dao*,  Chandrika  P.  Vyasarayani,  John  McPhee 

Department  of  Systems  Design  Engineering,  University  of  Waterloo,  200  University  Ave  West,  Waterloo,  Ontario,  Canada  N2L  3G1 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  1 1  July  201 1 
Received  in  revised  form 
13  September  2011 
Accepted  13  September  2011 
Available  online  21  September  2011 


Keywords: 

Galerkin 

Model  simplification 
Model  order  reduction 
Lithium-ion  battery 
Porous  electrode 
Hybrid  electric  vehicle 


In  this  paper,  effective  and  systematic  steps  in  the  mathematical  simplification  and  reduction  of  physics- 
based  lithium-ion  (Li-ion)  battery  models  to  improve  computational  efficiency  will  be  presented.  The 
battery  model  used  for  simulations  is  an  isothermal  model  proposed  by  Newman  and  Tiedemann  [1] 
and  Doyle  et  al.  [2]  which  incorporates  the  concentrated  solution  theory,  the  porous  electrode  theory, 
and  the  variations  in  electronic/ionic  conductivities  and  diffusivities.  The  simplified  model  is  formulated 
by  exploiting  the  nature  of  the  model  and  the  structure  of  the  governing  equations.  Simulations  show 
that  the  simplified  model  can  reduce  computational  time  significantly  while  still  retaining  the  accuracy 
compared  to  the  full-order  rigorous  model. 
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1.  Introduction 

Due  to  its  advantages  such  as  light  weight,  low  self-discharge 
rate,  and  high  specific  energy,  Li-ion  batteries  have  become  one  of 
the  most  popular  types  of  battery  in  various  applications  such  as 
portable  electronics  or  electric  vehicles.  The  charge-discharge  rate 
for  the  lithium-ion  battery  can  vary  from  1  to  2  C  (1  C  is  a  discharge 
rate  at  nominal  battery  capacity)  to  a  very  fast  pulse  discharge  up  to 
40-50  C  over  a  short  time  period  on  the  order  of  1 0-20  s  [3 ].  In  the 
automotive  field,  Li-ion  batteries  are  the  core  of  the  energy  source 
and  storage  in  new  plug-in  hybrid  electric  vehicles  (PHEVs)  as  well 
as  considered  in  many  second  generation  hybrid  electric  vehicles 
(HEVs).  The  performances  of  Li-ion  batteries  play  an  important  role 
in  vehicle  energy  management.  Hence,  battery  modeling  is  one  of 
the  most  important  tasks  for  electric  and  hybrid  vehicle  control. 
This  requires  a  model  that  can  simulate  in  real-time  in  order  to 
make  them  compatible  with  estimation  algorithms  embedded  in 
on-board  electronic  control  units.  For  example,  a  battery  manage¬ 
ment  system  in  an  HEV  has  to  estimate  the  battery  state  of  charge 
(SOC)  [4-6]  in  real-time  in  order  control  the  electrical  power  com¬ 
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ing  in  and  out  of  the  battery  as  well  as  to  prevent  the  battery  pack 
from  excessive  heating. 

A  full-order  battery  model,  however,  is  not  suitable  for  real¬ 
time  applications  as  it  usually  takes  hours  to  spatially  discretize 
the  system  via  a  finite  difference  method  and  solve  the  system 
numerically  as  a  collective  of  differential  equations  in  terms  of  the 
field  variables.  Therefore,  a  fast  and  reliable  approximate  model  is 
required.  For  automotive  applications,  a  simplified  battery  model 
has  to  be  carried  out  at  a  good  accuracy  while  ensuring  the  maxi¬ 
mum  computational  cost  reduction  to  achieve  an  efficient  system 
management.  Subramanian  et  al.  [7-9]  developed  a  real-time  sim¬ 
ulation  model  using  a  combination  of  perturbation  techniques, 
volume  averaging,  and  intuition-based  simplifications.  Although 
they  reported  that  the  computational  time  for  their  real-time  sim¬ 
ulation  model  for  a  single  process  was  around  100  ms,  to  derive 
the  lower-order  model  by  using  this  method  one  needs  to  carry 
out  preprocessing  and  have  a  priori  knowledge  of  the  behavior  of 
the  system  under  different  conditions,  which  makes  their  method 
less  flexible  than  desired.  Other  methods,  including  the  Chebyshev 
polynomial  method  [10,11],  residue  grouping  method  [12],  proper 
orthogonal  decomposition  method  [13],  and  Pade  approximation 
[14]  have  also  been  used  to  derive  reduced-order  models  for  Li-ion 
batteries. 

In  the  methods  using  Chebyshev  polynomials,  the  state  variables 
are  approximated  by  linear  combinations  of  several  Chebyshev 
polynomials,  and  then  an  approximate  model  is  projected  onto  a 
subspace  formed  by  these  orthogonal  Chebyshev  polynomials  to 
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Nomenclature 

ak 

specific  surface  area  of  electrode  k  (/< =p,  n)  (m2  m-3 ) 

bruggk 

Bruggman  coefficient  of  region  k  ( 1<=p ,  n) 

Ce.apprx 

assumed  solution  for  electrolyte-phase  concentra¬ 
tion  of  Li+  (molm-3) 

Ce,k 

electrolyte  concentration  in  region  k  (molm-3) 

Ce,k,  0 

initial  electrolyte  concentration  in  region  k 
(molm-3) 

Cs,k 

concentration  of  Li+  ions  in  the  intercalation  particle 
of  electrode  k  (molm-3) 

Cs,k,0 

initial  concentration  of  Li+  ions  in  the  intercalation 
particle  of  electrode  k  (molm-3) 

^s,k 

average  concentration  of  Li+  ions  in  the  intercalation 
particle  of  electrode  k  (molm-3) 

^s.k.surf 

concentration  of  Li+  ions  on  the  surface  of  intercala¬ 
tion  particle  of  electrode  k  (mol  m-3) 

D 

electrolyte  diffusion  coefficient  (m2  s) 

DsJ< 

Li+  ion  diffusion  coefficient  in  the  intercalation  par¬ 
ticle  of  electrode  k  (m2  s) 

F 

Faraday’s  constant  (C  mol-1 ) 

I 

applied  current  density  (Am-2) 

jk 

wall-flux  of  Li+  on  the  intercalation  particle  of  elec¬ 
trode  k  (molm-2  s) 

Kk 

intercalation/deintercalation  reaction-rate  constant 
of  electrode  k  (mol  mol-1  m3) 

L 

total  thickness  of  cathode-separator-anode  (m) 

h 

thickness  of  region  k  (m) 

n 

negative  electrode 

N 

number  of  node  points  for  Galerkin’s  approximation 

V 

positive  electrode 

Qs,k 

volume-averaged  concentration  flux  of  Li+ 
ions  in  the  intercalation  particle  of  electrode  k 
(molm-3  s-1) 

r 

radial  coordinate  (m) 

R 

universal  gas  constant 

Rc 

residual  function  for  concentration  of  Li+  in 
electrolyte-phase 

Rsj< 

radius  of  intercalation  of  electrode  k  (m) 

R<$> 

residual  function  for  electrical  potential  in 
electrolyte-phase 

s 

separator 

t+ 

Li+  transference  number  in  the  electrolyte 

T 

absolute  temperature  (I<) 

Uk 

open-circuit  potential  of  electrode  k  (V) 

X 

spatial  coordinate  (m) 

Greek  letters 

tk 

volume  fraction  of  region  k 

Gf,k 

volume  fraction  of  fillers  in  region  k 

m 

time-dependent  variable  of  i’th  basis  function  for 
electrolyte-phase  concentration  (s) 

0k 

dimensionless  concentration  of  Li+  ions  in  the  inter¬ 
calation  particle  of  electrode  k  (6k  =  cs  fc/cs  fc  max) 

K 

ionic  conductivity  of  electrolyte  (S  m-1 ) 

^eff J< 

effective  ionic  conductivity  of  the  electrolyte  in 
region  k  (S  m-1 ) 

Pi 

time-dependent  variable  of  i’th  basis  function  for 
electrolyte-phase  potential  (s) 

°k 

electronic  conductivity  of  solid  phase  of  electrode  k 
(Srrr1) 

^eff.k 

effective  electronic  conductivity  of  solid  phase  of 
electrode  k  (Sm-1) 

Oe  electrolyte-phase  potential  (V) 

^e.apprx  assumed  solution  for  electrolyte-phase  potential  (V) 
Os  solid-phase  potential  (V) 


form  a  reduced-order  model,  which  can  then  be  solved  for  the 
unknown  coefficients  in  the  truncated  expressions.  Smith  et  al. 
[15]  developed  a  control-oriented  one-dimensional  (ID)  electro¬ 
chemical  model  by  using  the  method  of  residue  grouping.  Their 
transfer  functions  are  represented  by  a  truncated  series  of  grouped 
residues  with  similar  eigenvalues.  Cai  and  White  [13]  proposed  an 
approach  based  on  proper  orthogonal  decomposition  for  tackling 
the  problem  by  using  two  steps  of  approximation:  partial  differen¬ 
tial  equation  (PDE)  discretization  and  truncation  of  the  number  of 
orthogonal  modes.  Cai  and  White  showed  that  the  order-reduced 
model  simulated  seven  times  faster  than  the  full-order  model  for 
a  similar  level  of  accuracy.  The  electrode-averaging  method  was 
used  by  Speltino  et  al.  [16]  who  made  several  simplifications  such 
as  neglecting  solid  concentration  distribution  and  assuming  con¬ 
stant  electrolyte  concentration.  As  a  result,  their  model  simulated 
fast,  but  with  a  heavy  loss  of  information. 

In  this  paper,  an  efficient  method  for  reducing  the  order  of 
Li-ion  battery  models  using  LiCo02  and  LiC6  electrodes  derived 
from  the  porous  electrode  theory  will  be  discussed.  The  simplified 
model  in  this  paper  uses  the  volume-average  integration  proposed 
in  Subramanian  et  al.  [7],  Wang  et  al.  [17],  and  Gu  et  al.  [18]  for 
approximating  the  Li+  concentration  in  the  solid  phase  in  the  elec¬ 
trode  material.  For  modeling  the  Li+  concentration  and  electrical 
potential  in  the  electrolyte  phase,  Galerkin’s  approximation  will  be 
used  under  the  assumption  of  a  galvanostatic  discharge  condition. 
Maple  14  was  used  to  generate  the  mathematical  model  and  per¬ 
form  many  of  the  model  reductions  and  simplifications.  The  paper 
begins  with  a  brief  overview  of  the  intercalated  Li-ion  models  and 
its  mathematical  governing  equations,  followed  by  a  discussion  on 
the  use  of  the  volume-average  and  Galerkin’s  method  to  simplify 
the  model.  Simulation  results  and  a  comparison  between  the  sim¬ 
plified  and  full-order  models  will  also  be  provided  to  demonstrate 
the  effectiveness  of  the  proposed  battery  model  reduction  method. 

2.  Lithium-ion  battery  model  based  on  porous-electrode 
theory 

2.1.  Overview 

Most  of  the  current  rigorous  Li-ion  battery  models  are  derived 
from  the  porous  electrode  and  concentrated  solution  theories  pro¬ 
posed  by  Newman  and  Tiedemann  [1]  and  Doyle  et  al.  [2]  which 
mathematically  describe  charge/discharge  and  species  transport  in 
the  solid  and  electrolyte  phases  across  a  simplified  ID  spatial  cell 
structure.  This  ID  model  of  a  Li-ion  battery  considers  dynamics 
along  only  one  axis  (the  horizontal  x-axis)  and  neglects  the  dynam¬ 
ics  along  the  remaining  two  axes  (y-axis  and  z-axis)  [1,2,19-23]. 
This  approximation  is  applicable  to  most  cell  structures  as  the 
length  scale  of  a  typical  Li-ion  cell  along  the  x-axis  is  on  the  order 
of  100  |jim,  whereas  the  length  scale  for  the  remaining  two  axes  is 
on  the  order  of  100,000  p,m  or  more  [24]. 

There  are  four  main  components  in  a  typical  Li-ion  cell  as  shown 
in  Fig.  1 :  the  negative  composite  electrode  connected  to  the  nega¬ 
tive  terminal  of  the  cell,  the  positive  composite  electrode  connected 
to  the  positive  terminal  of  the  cell,  the  separator,  and  the  elec¬ 
trolyte.  The  negative  electrode  contains  lithium  stored  in  the  lattice 
sites  made  from  graphite,  usually  in  the  form  of  LixC6.  The  posi¬ 
tive  electrode  can  have  various  chemistries,  usually  a  metal  oxide 
or  an  inter-metallic  oxide  such  as  LixMn204  or  LiyCo02.  During 
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Fig.  1.  Anatomy  of  Li-ion  cell. 


discharge,  Li+  ions  diffuse  to  the  surface  of  LixC6  active  material  par¬ 
ticles  (solid  phase)  in  the  negative  electrode  where  they  undergo 
electrochemical  reactions  and  transfer  into  a  liquid  electrolyte  solu¬ 
tion  (solution  phase).  The  positively  charged  Li+  ions  travel  through 
the  electrolyte  solution  via  diffusion  and  ionic  conduction  to  the 
positive  electrode  where  they  react  and  diffuse  towards  the  inner 
regions  of  metal  oxide  active  material  particles  (solid  phase).  The 
process  is  similar  for  charging,  except  that  the  Li+  ions  move  in  the 
opposite  direction  from  the  lattice  sites  in  the  positive  electrode 
and  enter  the  lattice  sites  in  the  negative  electrode.  This  process  is 
called  intercalation  and  this  is  the  reason  why  the  Li-ion  battery  is 
usually  referred  to  as  the  rocking  chair  model.  The  porous  separator 
serves  as  an  electric  insulator  which  forces  electrons  to  follow  an 
opposite  path  to  the  ions  through  an  external  circuit  or  load.  The 
separator,  however,  allows  the  Li+  ions  to  pass  through  it  during 
battery  operation. 


2.2.  Transport  in  solid  phase 

In  this  section  we  provide  a  brief  summary  of  the  governing 
equations  for  a  lithium-ion  battery  model  derived  from  the  porous 
electrode  theory  [1,2,19].  According  to  this  theory,  the  lithium  is 
considered  as  existing  in  two  disjoint  states  called  phases:  the  solid 
phase  in  the  electrode  material  and  the  liquid  phase  in  the  dissolved 
state  in  the  electrolyte.  In  this  model,  the  solid  and  electrolyte 
phases  are  treated  as  superimposed  continua  without  regard  to 
microstructure.  According  to  Fields  laws  of  diffusion,  the  solid- 
phase  Li+  concentration  in  a  single  spherical  active  material  particle 
can  be  described  as 

9cSjk(x,  r,  t)  Ds  k  9  /  2  dcsk(x,  r,  t)\ 

9t  r2  dr  {  dr  J  1  j 

with  boundary  and  initial  conditions 


where  /<  =  p  for  the  positive  electrode  and  /<  =  n  for  the  negative 
electrode.  The  definitions  of  other  symbols  can  be  found  in  the 
Nomenclature  section. 


2.3.  Transport  in  electrolyte 


The  Li+  concentration  in  the  electrolyte  phase  changes  due  to 
the  changes  in  the  gradient  diffusive  flow  of  Li+  ions.  Therefore,  it 
can  be  shown  that 


dce,k(x,  t ) 

at 


d_ 

dx 


3ce,k(x,  t)\ 

dx  ) 


+  -t+)Mx,  t) 


(5) 


In  this  equation,  k  =  p  for  the  positive  electrode,  k  =  s  for  the  sepa¬ 
rator,  and  k  =  n  for  the  negative  electrode. 

The  boundary  conditions  for  Eq.(5)  capture  the  fact  that  the 
fluxes  of  the  ions  are  zero  for  all  time  at  the  current  collector.  Since 
the  flux  is  proportional  to  the  concentration  gradient  at  the  positive 
and  negative  current  collectors,  we  have 

n  9ce,p  n  9ce,n  ,  n  ,c\ 

— Z^eff,p  g  x=®  —  — ^eff,n  I  x=L  —  0  (6) 

We  also  need  four  additional  boundary  conditions  at  the 
electrode-separator  interface.  These  boundary  conditions  are 
obtained  from  continuity  of  the  flux  and  concentration  of  the  elec¬ 
trolyte  at  the  electrode-separator  interface  as 


n  ^Ce,p  1  n  ^Ce’s  1 

eff’p  dx  Ix=1p"  “  eff's  dx  lx=i? 

(7) 

n  ^Ce’s  1  n  ^Ce’n  1 

°eff,s  dx  x={Lp  +Ls)~  ~  ^eff’n  dx  'X=(b+Ls)+ 

(8) 

ce,p  1  x=L~  =  ce,slx=L+ 

0) 

Ce,slx=(Lp+Lsr  =  Ce,nlx=(Lp+Ls)+ 

(10) 

Eq.  (5)  must  also  satisfy  the  initial  condition 


ce,k(x,  0)  =  ceX0  (11) 

The  effective  diffusion  coefficient  Deff  fc  is  calculated  from  a  refer¬ 
ence  coefficient  using  the  Bruggman  relation  Deff  k  =  De£rugg/<  that 
accounts  for  the  tortuous  path  that  Li+  ions  follow  through  the 
porous  media.  In  this  expression,  D  is  the  electrolyte  diffusion  coef¬ 
ficient  which  varies  with  electrolyte  concentration.  However,  D  will 
be  approximated  as  constant  within  this  paper  for  simplification. 
The  specific  electrode  surface  area  can  be  expressed  in  terms  of  the 
porosity  of  the  electrode  as 

ak  =  n —  (1  -  ek  -  €f,k)  02) 

,k 


2.4.  Electrical  potentials 


Charge  conservation  in  the  solid  phase  of  each  electrode  can  be 
described  by  Ohm’s  law 

aeff,k - ^  =  akFJk(x >  0  03) 

with  boundary  conditions  at  the  current  collectors  being  propor¬ 
tional  to  applied  current  density 


9cs  u 

—Ds,k  lr=0  =  0 

~Ds,k  lr=fts =Jk(xi  0 

Cs,k(x,  r’  °)  =  cs,k,0 


(2) 

(3) 

(4) 


90s,p  90s>n  T 

^  lx=0  —  —  ^eff.n  ^3  lx=L  —  1 


^eff,  |: 


90 


s,p 


dx 


lx=Lp  —  ^eff,!! 


dx 

a^s,n 

dx 


lx— Lp+Lg  ^ 


(14) 

(15) 


where  the  current  density  /  is  related  to  the  applied  current  i  and 
the  surface  area  A  of  the  electrode  as  /  =  i/A. 
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The  effective  electronic  conductivity  can  be  expressed  in  terms 
of  the  porosity  of  the  electrode  as 

aei f,k  =  °k  (l  -  €k  -  €f ,k)  06) 

Combining  Kirchhoff  s  law  with  Ohm’s  law  in  the  electrolyte 
phase  yields 


30s?k(x,  t)  3 ®e,k(x>  0 

■ffeff’k  ax  “*eff’k  ax 

2icefUk(x,t)RT  .3lnce,k 

H - F - 1 1  ~  l+) - q -  — 


(17) 


with  9k  being  given  by 

ek(x,  t)  =  Cs-k’surf(x’ 

^s,k,max 


(25) 


The  Butler-Volmer  equation  describing  the  relationship 
between  the  current  density,  concentrations,  and  over-potential 
is  given  by  Newman  and  Thomas-Aleya  [21  ] 


Jkfai  0  —  ^ki^s,k, max 

'0.5  F 


exp  | 


RT 


Ms, 


-  Cs,fc,surf)°'5(Cs,k,surf)0'5Cej(X,  t) 

k(x,  t))  -  exp  (-^*(*.  t)) 


(26) 


Both  t+  and  /ceff tk  are  usually  functions  of  electrolyte  concen¬ 
tration,  but  t+  is  usually  approximated  as  a  constant.  The  effective 
diffusion  conductivity  can  be  calculated  from  the  Bruggman  rela¬ 
tion.  In  this  paper,  the  conductivity  for  the  liquid/salt/polymer 
system,  consisting  of  a  2:1  volume  mixture  of  ethylene  carbonate 
(EC)  and  dimethyl  carbonate  (DMC),  curve-fitted  by  Doyle  et  al.  [25] 
will  be  used  in  the  following  form 


In  summary,  the  equations  that  need  to  be  solved  are  (1),  (5), 
(13),  (17),  (23),  (24),  and  (26).  The  battery  model  is,  therefore,  a 
mixed  system  with  14  nonlinear  partial  differential  algebraic  equa¬ 
tions  (PDAEs)  with  14  unknowns,  which  are:  cS)P,  cs,n,  ce,p,  ce,s,  ce,n, 

3>s,p,  3>s,n,  Oe.p,  Oe.s,  Oe,n,  UP>  U*  Jp»  ^djn. 

3.  Model  simplification 


/ceff >k(x,  t)  =  e£ruggk[4.1253  x  10  2  +  5.007  x  10  4cek{x ,  t )  3.1.  Solid  phase  concentration 


-4.7212  x  10_7c2  k(x,  t)  +  1.5094  x  10~10c^k(x,  t) 

-  1.6018  x  10 ~14c4k{x,t)]  (18) 


Since  we  can  only  measure  potential  differences,  the  boundary 
conditions  of  Oe/<  is  arbitrary.  We  set  Oe,p(0+,  t)  =  0  at  the  posi¬ 
tive  electrode  current  collector  interface.  The  remaining  boundary 
conditions  follow  from  continuity  of  Oe  k  that 


30, 


^eff,P  ^  l*=0  —  ^eff, n 


3Qe,n 

dx 


I  x=L  =  0 


(19) 


^eff,P 


9^e,p, 

dx  Ix=Lp" 


^eff,  s 


3<*>e,S 

dx 


lx=i+ 


(20) 


For  the  solid-phase  concentration,  we  are  only  interested  in 
the  solid-phase  Li+  concentration  on  the  surfaces  of  the  lithium 
particles  for  electrochemical  behaviors  and  average  concentration 
for  state  of  charge  calculation.  Therefore,  it  is  desirable  to  only 
extract  the  equations  for  these  two  quantities  without  having  to 
solve  the  entire  two  PDEs  for  the  solid-phase  concentration.  In 
the  process,  the  dependencies  on  the  radial  dimension  r  are  elimi¬ 
nated. 

Simplified,  yet  accurate,  equations  for  surface  and  average  con¬ 
centrations  can  be  obtained  using  the  polynomial  approximation 
and  volume-average  integration  suggested  by  Subramanian  et  al. 
[7].  The  procedures  can  be  started  by  choosing  a  concentration  pro¬ 
file  inside  a  particle  at  a  known  x  position.  In  this  case,  we  choose 
a  three-variable  model 


30e,s  30e,i 

ICeff,s ~dx~ Ix=(lp+Ls)_  “  ~Kef^n~dx~ 


=(Lp  +Ls)+ 


(21) 


It  should  be  noted  that  the  solid-phase  potential  Oe  fe(x,  t)  and 
current  density  Jk{x ,  t)  do  not  exist  in  the  separator  region  and, 
therefore,  the  terms  containing  these  variables  are  eliminated  from 
Eqs.  (5)  and  (17)  for  the  separator  region. 


2.5.  Butler-Volmer  kinetics 


cs,k(r ,  t )  =  a(t)  +  b(t)-=-  +  d(t)—r— 


Rlk 


Kk 


(27) 


Substituting  the  above  polynomial  function  into  Eq.  (1 )  yields 


da(t)  r 
~dF~  + 


2  db{t)  t 


dd(t ) 


Rlk 


dt 


R4s.k 


dt 


-2t+*t)+,V(+° 

(28) 


The  molar  flux  Jk(x ,  t)  depends  on  the  concentration  cs  k  of  Li+  in 
the  electrode  k,  the  concentration  of  Li+  in  the  electrolyte,  and  the 
intercalation  over-potential  /xs/<(x,  0  through  the  Butler-Volmer 
equation  [21  ].  This  over-potential  can  be  described  as 


M s,k(x,  t)  =  <Ds,k(x.  t)  -  <J>e,k(x,  t)  -  Uk(9k(x,  t))  (22) 

where  the  open-circuit  potential  equations  for  positive  (LiCo02) 
and  negative  (LiC6)  electrodes  were  curve-fitted  from  experimental 
data  [2]  and  have  the  form 


This  equation  automatically  satisfies  the  boundary  condition  at 
r  =  0,  therefore  we  only  need  to  make  sure  the  equation  also  satisfies 
the  second  boundary  condition  at  r=Rsk.  This  corresponds  to 


2Ds,km  4Ps>fcd(t) 
Rsk  Rs,k 


(29) 


The  three  unknowns  a(t),  b(t ),  and  d(t)  can  be  solved  in  terms  of 
the  average  concentration  cs  k(t),  surface  concentration  cSik, surfM- 
and  an  auxiliary  term  called  volume-averaged  concentration  flux 


Up(x,  t)  = 


-4.656  +  88.66902  _  401.11904  +  342.90906  -  462.47106+433.43401° 
-1.0  +  18.93302  -  79.53204  +  37.31106  -  73.0830®  +  95.960’° 


(23) 


Unix,  t )  =  0.7222  +  O.13870n  +  O.O290°  5 


0.0172  0.0019 

~^r  +  ^p~ 


qStk(t)  as  (detailed  solution  can  be  found  in  Subramanian  et  al.  [7]) 

39  _  35_ 

a(0  —  ^cs,k,surf(0 —  2qsk{t)  ~^Csk{t) 

b(t)  =  —  35cs  ?fc?surf(t)  +  1  Oqsk(t)  +  35cs  ,k(t) 
d(t)  =  -^-cs  k  surf(t)  -  7qs  k(t)  -  -^-cs  k(t) 


+  0.2808  exp(0.90  -  150„)  -  0.7984 
x  exp(O.44650„  -  0.4108) 


(24) 


(30) 
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The  solid-phase  concentration  now  can  be  expressed  in  terms 
of  the  surface  concentration  and  average  concentration  by  substi¬ 
tuting  (30)  into  (28) 


39  _  35_ 

Cs,k(ri  0  —  “l“^s,k,surf(0  —  Bqsk(t)Rsk  k(t)[  — 35cs ?fc?surf(t) 


T  1 0^s,/c(f )(f )^s,k  +  35cs  /<(£)]  - 

s,k 


105 


105_ 


/I  ^s,/c,surf(0  7qsk{t)Rs,k  a  Cs,k(t) 


,  -=-  (31) 


By  applying  the  volume-average  integration  for  the  original  PDE 
and  its  differentiation,  we  can  obtain  the  two  ordinary  differential 
equations  (ODEs) 


d-cs  k(t)  +  =  o 

^s,fe 


D, 


+  30 -  0 


45Jk(t) 


(32) 

(33) 


The  third  equation  can  be  obtained  by  evaluating  the  boundary 
condition  at  r=Rsk.  This  gives 


time  rji(t).  It  can  be  easily  verified  that  the  assumed  solution  func¬ 
tion  satisfies  the  boundary  conditions  as  well.  The  ce,o  is  included  to 
ensure  the  assumed  solution  satisfies  the  initial  condition  given  in 
Eq.  (11).  The  first  two  terms  in  the  assumed  solution  in  (37)  give  a 
straight  line  function  of  x  that  interpolates  the  boundary  values.  The 
terms  in  the  summation  contribute  nonlinearities  to  the  solution. 

Substituting  the  assumed  solution  into  the  PDE  in  (35)  gives 


N 

RC(X,  t)  =  €^2  («iM 
i= 1 


dt 


d2  cq(x) 


+  a(l-t+l/«0 


(38) 


This  function  is  known  as  the  residual  In  the  Galerkin  method, 
we  replace  the  condition  that  the  residual  should  be  approximately 
zero  by  the  condition  that  the  residual  should  be  orthogonal  to  the 
set  of  basis  functions.  That  is,  for  i  =  1 , . . .,  N  we  multiply  the  residual 
by  the  basis  function  cos  (ittx/L)  and  integrate  over  x,  and  set  the 
result  to  zero.  For  convenience,  we  define  the  following  operator: 

Definition  1.  Inner  product  Consider  the  real  functions  whose 
domain  is  the  closed  interval  [a,  b].  We  define  the  inner  product  of 
two  functions /(x,  t)  and  g(x)  as  follows 


35^[cs,k,surf(t)  -  cs>k(t)]  -  8DS  kqs  k(t)  =  -Jk(t)  (34) 

Ks,k 

Solving  the  ODEs  in  (32),  (33),  and  (34)  simultaneously  using 
a  numerical  solver  gives  the  surface  and  average  concentration 
profiles  for  the  Li+  concentration  in  solid  phase. 

The  next  sections  will  present  the  main  contribution  of  this 
paper,  the  use  of  Galerkin’s  method  and  analytical  techniques  to 
reduce  the  PDEs  describing  the  electrolyte  phase  concentration 
(cetfc)  and  electrical  potential  equations  (Os <k  and  Oe  fe)  to  ODEs. 

3.2.  Electrolyte  phase  concentration 


( f,g)a=  f  f(x,t)g{x)dx  (39) 

J  a 

Using  this  inner  product  operator,  we  can  write  the  integration 
as  follows 

(Rc,  a,)^  +  (Rc,  ai)[p+Ls  +  (Rc,  a,)|p+ls  =0,  i  =  1 , . . . ,  N  (40) 

In  the  above  equation,  we  substitute  (e  =  en,  Deff  =  Deffn,  a  =  an) 
into  the  first  term,  (e  =  es,  Deff  =  A;ff,s»  a  =  as )  into  the  second  term, 
and  (e  =  ep,  Deff  =  DeffiP,  a  =  ap )  into  the  last  term.  In  this  way  we 
obtain  a  set  of  N  linear  ODEs  that  only  contains  the  time-dependent 
functions  rji(t).  These  ODEs  can  be  expressed  using  the  matrix  form 


The  Li+  concentration  in  the  electrolyte  phase  can  be  approx¬ 
imated  by  applying  Galerkin’s  approximation  to  Eq.(5).  The  first 
step  of  the  Galerkin  method  is  to  choose  a  basis  function  that  satis¬ 
fies  all  the  boundary  conditions  (6)-(7),  and  (8).  For  simplicity,  we 
approximate  the  three  PDEs  across  three  regions  (i.e.,  anode,  sepa¬ 
rator,  and  cathode  regions)  by  one  single  PDE  with  the  x-dimension 
spanning  from  0  to  L 


/  m(t)\ 

(  »7i(t)\ 

Ut)  =-I^l 

r]2(t) 

\  >1nU)  ) 

where  [Ac]  and  [Bc] 

are  NxN 

+  [Bc] 


are  obtained  during  the  integration  process  in  (40). 


(41) 


3 Ce(x,  t) 

3t 


d_ 

3x 


3 ce(x,  t) 
3x 


+  a(l  -  t+)J 


(35) 


The  only  differences  in  the  three  parts  of  this  PDE  are  the  phys¬ 
ical  parameters  in  the  three  regions  and  the  boundary  conditions 
at  the  cathode-separator  and  separator-anode  interfaces.  This  can 
be  taken  into  account  when  applying  the  least-square  integration 
and  will  be  shown  in  the  next  steps. 

Using  the  criteria  discussed  above,  we  can  choose  the  sinusoidal 
basis  function 

«;(*)  =  c< °s(^)  (36) 

We  can  then  define  the  assumed  solution  as  the  sum  of  the  basis 
functions  in  the  form 


3.3.  Solid  phase  potential 

For  simplification,  we  assume  that  the  current  density  Jk  is  con¬ 
stant.  Based  on  this  assumption,  we  can  solve  for  the  closed-form 
solution  for  the  solid-phase  potential.  Substituting  Jp  =  //apFLp  into 
Eq.  ( 1 3 )  and  integrating  Eq.  ( 1 3 )  twice  over  x  gives  a  quadratic  equa¬ 
tion  for  the  positive  solid-phase  potential 

4>s,p(x,  t)  =  — laPFJp*  +gp(t)x+/p(t)  (42) 

1  ^eff,  p 

where  gp(t)  and /p(t)  are  the  quantities  produced  by  the  indefinite 
integration.  The  expression  for  gp(t)  can  be  determined  by  evalu¬ 
ating  Eq.(42)  at  the  boundary  condition  at  x  =  0  shown  in  Eq.(14). 
This  gives 


N  N 

Ce,apprx(X,  t)  =  Ce,0  +  ^[a, •(*)>?, -(f)]  =  Ce, 0  +  [cOS  (^)  rjj(t) 

1  =  1  1  =  1 

(37) 

where  N  is  the  number  of  terms  or  node  points.  Each  term  of  the  sum 
is  the  product  of  a  given  basis  function  and  an  unknown  function  of 


gP(t) = 


i 

^eff,  p 


which,  when  substituted  into  Eq.(42)  yields 


4>s,p(x,  t)  = 


1  Qpf/p*2 
3  fJetr.p 


I 

^eff,  p 


X+/p(t) 


(43) 


(44) 
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Similarly,  by  integrating  Eq.  (13)  twice  over  x  and  evaluating  its 
result  at  the  boundary  condition  at  x  =  L  the  closed  form  solution  for 
the  solid-phase  potential  on  the  negative  electrode  can  be  obtained 
as 


<*>s,n(x,  t)  = 


_  \  ClnFJnx2  +  anf/n(^s  +  ^p)  ^  +  j 
2  cre  ff,n  ~  “ 


(45) 


From  Eqs.  (44)  and  (45),  instead  of  solving  for  Os,p(x,  t)  and 
Os,n(x,  t),  we  can  now  substitute  these  two  quantities  into  the 
electrolyte  phase  potential  (discussed  in  the  next  section)  and 
Butler-Volmer  Eq.  (26)  and,  at  a  known  x  position,  we  can  solve 
for  the  two  unknowns /p(t)  and  /n(t). 


3.4.  Electrolyte  phase  potential 


Similar  to  the  process  used  for  the  electrolyte-phase  concen¬ 
tration,  we  can  apply  Galerkin  method  to  the  electrolyte  phase 
potential  to  obtain  approximate  ODEs  for  this  PDE.  In  a  similar 
manner,  we  use  a  single  PDE  to  represent  the  three  PDEs  in  three 
different  regions  as  follows 

30s(x,  t)  30e(x,  t)  2/ceff(x,  t)RT 
~CTeff  dx  “ Ke<<  dx  + - F - 


x  (1  -  4) 


3lnce(x,  t) 
dx 


=  1 


(46) 


Fig.  2.  Open-circuit  voltage  for  positive  electrode. 


rji(t )  which  come  from  the  integration  of  x(x,  t)  and  In  ce(x,  t)  while 
the  last  two  matrices  are  constant. 


In  a  similar  manner,  we  can  choose  the  basis  function 
A (x)  =  cos  (i'ttx/L)  and  write  the  approximate  function  as 


4>e,apprx(X,  0  =  ^[AMaWI  =  ^  [cOS  ^  — j  A(0 

i=l  i=l 

and  then  the  residual  by  substituting  Eq.  (47)  into  Eq.  (46) 

N 


K0(X,  t)  =  -CTeff 


3 <&s(x,t)  ,  ,,  s^fdPiix) 


dx 


‘^eff 


2/ceff(x,  t)f?T 


^  V  dx 

i=i  x 

^  3lnce(x,  t) 
dx 


Pi(t) 


(47) 


(48) 


The  set  of  N  equations  can  be  obtained  by  the  following  integra¬ 
tion 


(R<fi,  Pi)(F  +  (fy»  Pi)LLPp+Ls  +  Pi'lLp+Ls  =  °>  «  =  (49) 

In  the  above  equation,  we  substitute  appropriate  parameters 
for  each  region  and  use  Eq.  (44)  for  the  positive  electrode  solid- 
phase  potential  and  Eq.(45)  for  the  negative  electrode  solid-phase 
potential.  The  electrolyte-phase  concentration  ce(x,  t)  inside  the 
ln(.)  function  can  also  be  replaced  by  its  approximate  expression 
in  Eq.  (37)  and  evaluated  at  the  nodal  x  positions.  Doing  so  we  obtain 
a  set  of  N  nonlinear  algebraic  equations  (there  are  no  ODEs  because 
there  are  no  time  derivatives  in  the  original  PDE)  containing  a  mix¬ 
ture  of  the  time-dependent  functions  a(0,  A‘(f)>/p(f)>  and/n(t)  in 
the  following  matrix  form 


[\{rji{t))] 


(  Pi(0\ 

A>(0 


+ 


\  Av(0  / 


/  mix) 

mi*) 

\mit) 


+  [D(j)\ 


//p(0n 

/p(0 

\fn(t)y 


+  [F(f)] 


(50) 


in  which  the  inverse  of  rji(t)  appears  due  to  the  differentiation  of  the 
logarithm  terms.  The  first  two  square  matrices  are  the  functions  of 


3.5.  Open-circuit  voltage  equations 


The  order  of  the  open-circuit  voltage  equations  for  the  two 
electrodes  given  in  (23)  and  (24)  can  be  reduced  to  increase  com¬ 
putational  speed  by  eliminating  the  high-order  terms  using  a 
nonlinear  curve-fitting  technique.  After  trying  different  function 
combinations,  we  have  come  up  with  simpler  forms  for  these  func¬ 
tions.  The  original  open-circuit  equations  can  be  replaced  by  the 
following  simple  expressions 


-4.875  +  5.8393p  -  1.5O70jj  +  0.531^ 
Up(t)= -  p  p  p 


l/n(t)  =  0.15-0.106>n  + 


Op  -  1.005 
0.00778 


On 


(51) 

(52) 


For  comparison,  the  open-circuit  voltages  in  Eqs.  (51)  and 
(52)  plotted  along  side  with  the  original  equations  are  shown  in 
Figs.  2  and  3. 


3.6.  Solving  battery  equations 


Using  the  procedures  discussed  above,  we  obtain  a  set  of  2N  + 1 0 
differential  algebraic  equations  (DAEs)  that  consists  of  the  depen- 
dent  variables  ^-(t),  p,(t),  fp(t),  /n(t),  Up(t),  Un(t),  cSik(t),  qStk(t), 
cs,k,surf(f)  and  two  independent  variables  x  and  t.  These  vari¬ 
ables  correspond  to  the  differential  and  algebraic  expressions  in 
(32)-(34),  (41),  (50)-(52),  and  one  (26)  for  each  of  Jp  and  Jn  with 
the  assumption  of  uniform  reaction  rate. 

We  are  mainly  interested  in  calculating  the  battery  voltage 
Veen (t)  given  by  the  following  relation 

^cell(t)  =  4>s,p(0,t)-Os,n(L,t)  (53) 

This  equation  results  in  two  sets  of  DAEs  evaluated  at  x  =  0  and 
x  =  L.  These  DAEs  are  only  functions  of  time  and  can  be  solved 
numerically  if  the  initial  conditions  are  known. 

The  initial  conditions  for  all  variables  are  needed  to  solve  the 
system.  For  the  solid-phase  concentration  variables  cs  k{t)  and 
cs,k,surf(f)  the  initial  conditions  are  known  and  are  equal  to  cs<k0.  The 
initial  value  for  the  third  variable  qS)k(t)  in  the  solid-phase  concen¬ 
tration  equations  can  be  determined  by  evaluating  Eqs.  (32)-(34) 
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_ n _ 

□  Original - Curve-fitted 


Fig.  3.  Open-circuit  voltage  for  negative  electrode. 


4-  Simplified  model:  1C  discharge  □  Simplified  model:  0.5C  discharge 
—  ■  —  Rigorous  model:  1C  discharge - Rigorous  model:  0.5C  discharge 


Fig.  4.  Discharge  curves  for  1  C  (30  A  rrr 2 )  and  0.5  C  rates  at  N  =  4. 


Fig.  5.  1  C  discharge  voltage  curve  comparison  between  the  rigorous  model  and  the 
simplified  model  at  different  number  of  node  points. 


Fig.  6.  Voltage  curves  for  different  discharge  rates  at  N  =  4. 


at  cs>fe(t)  =  cs  k  surf{t)  =  cs  k  0.  In  the  DAEs  obtained  from  Galerkin’s 
method  in  Eqs.  (41 )  and  (50),  we  can  solve  for  the  initial  conditions 
for  r]i(t)  and  pf(t)  by  first  expanding  the  initial  conditions  of  the 
assumed  solution  in  terms  of  a  series  solution  as 

N 

Wapprx(*?  °)=£  Ti(xMO)  (54) 

i=i 

where  T(x)  represents  the  x-dependent  part  and  f  (0)  represents  the 
initial  value  of  the  time-dependent  part  of  Eqs.  (41)  and  (50).  We 
can  then  obtain  the  values  for  §(0)  (i.e.,  ^(0)  or  a(0))  by  multiplying 
the  above  equation  with  T(x)  and  solving  the  resulting  equations 

N 

<wapprx(x,  o).  r,(x)>&  -  ]r^(0)<r»,  r,(x)>£  =  o,  i  =  i , . . . ,  n 

J=1 

(55) 


Fig.  7.  Electrolyte-phase  concentration  of  Li+  at  the  current  collector  and  electrode¬ 
separator  interfaces  at  1  C  discharge  rate. 
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Fig.  8.  Electrolyte-phase  concentration  of  Li+  at  the  current  collector  and  electrode¬ 
separator  interfaces  at  1  C  discharge  rate. 

The  initial  values  for  the  remaining  variables  Up(t ), 

and  Un(t)  can  be  obtained  by  solving  all  the  algebraic  expressions 
with  the  known  initial  conditions. 

3.7.  Simulation  results 

Discharge  behavior  is  the  primary  goal  of  the  simulation  per¬ 
formed  in  this  paper.  The  parameters  used  in  the  simulations  are 
obtained  from  Subramanian  et  al.  [9]  and  are  listed  in  Table  A.l. 
Fig.  4  shows  complete  discharge  voltage  curves  at  1  C  (30 Am-2) 
and  0.5  C  rates  of  galvanostatic  discharge  with  the  number  of  node 
points  being  N  =  4.  Fig.  5  compares  the  voltage  curves  resulting 
from  the  simplified  model  at  different  numbers  of  node  points 
and  the  full-order  rigorous  model  solved  using  a  finite  difference 
approach  at  1  C  discharge  rate.  The  plots  for  the  rigorous  model 
were  extracted  from  the  work  done  by  Subramanian  et  al.  [9].  It 
can  be  seen  that  even  with  the  number  of  node  points  N=  2,  there 
is  a  good  agreement  between  the  simplified  model  and  the  rigorous 
model  and  this  was  improved  further  at  N  =  4.  The  number  of  DAEs 
that  are  solved  simultaneously  using  Galerkin’s  approach  is  2N  + 10 
while  the  rigorous  model  needs  to  be  discretized  into  hundreds  of 


Fig.  9.  Electrolyte-phase  potential  at  the  current  collector  and  electrode-separator 
interfaces  at  1  C  discharge  rate. 


equations  to  obtain  a  good  accuracy.  For  example,  the  results  for  the 
full-order  model  shown  in  Figs.  4  and  5  were  obtained  by  discretiz¬ 
ing  the  PDEs  governing  the  battery  model  into  600  DAEs  shown 
in  the  work  of  [9].  The  simplification  also  results  in  a  great  reduc¬ 
tion  in  simulation  time  for  the  approximate  model.  Specifically,  it 
took  94  ms  on  average  to  simulate  the  system  at  N  =  4  in  Maplesoft’s 
flagship  software,  Maple  14,  on  a  Dell™OptiPlex  2.9GFIz  desktop 
computer. 

The  battery  potentials  for  low  and  high  discharge  rates  are  plot¬ 
ted  in  Fig.  6.  Since  no  results  for  high-rate  discharge  that  can  be 
used  for  comparison  were  reported  in  the  work  by  Subramanian 
et  al.  [9],  the  consistency  between  the  results  over  a  range  of  dis¬ 
charge  rates  shown  in  Fig.  6  indicates  that  the  simplified  model 
works  as  expected. 

The  simplified  model  is  also  able  to  predict  the  other  electro¬ 
chemical  variables  such  as  the  concentration  of  Li+  in  the  electrolyte 
(ce)  as  shown  in  Figs.  7  and  8,  and  the  electrolyte  phase  poten¬ 
tial  (<I>e)  as  shown  in  Fig.  9.  The  electrolyte-phase  concentration 
and  potential  determined  by  the  simplified  model  are  plotted  at 
t=  10,  t= 20,  t=  100,  and  t=1000.  It  can  be  seen  that  there  is  no 
discontinuity  in  the  slopes  of  the  curves  at  the  positive  electrode¬ 
separator  and  separator-negative  electrode  boundaries  because  the 


Table  A.1 

Battery  parameters. 


Symbol 

Unit 

Positive  electrode 

Separator 

Negative  electrode 

Sirr1 

100 

100 

€f ,k 

0.025 

0.0326 

€k 

0.385 

0.724 

0.485 

brugg 

4 

4 

4 

Ds,k 

m2  s_1 

1  x  10-14 

3.9  xlO”14 

D 

m2  s_1 

7.5  x  10-10 

7.5  xlO-10 

7.5  x  10-10 

I<k 

mol  (mol  rrr3)-1-5 

2.344  x  10-11 

5.0307  x  10-11 

Cs,k,max 

mol  m3 

51,554 

30,555 

Cs,k,0 

mol  m3 

0.4955x51,554 

0.8551  x  30,  555 

Ce.O 

mol  m3 

1000 

1000 

1000 

Rs,k 

m 

2.0  xlO-6 

2.0  x  10“6 

Lk 

m 

ID 

o 

X 

o 

00 

25  x 10“6 

00 

00 

X 

o 

O) 

P-SEl 

Q  m2 

0 

u 

0.363 

0.363 

0.363 

F 

Cmol-1 

96,487 

R 

J  mol-1  K_1 

8.314 

T.-S.  Dao  et  al.  /  Journal  of  Power  Sources  198  (2012)  329-337 


337 


three  PDEs  in  three  regions  are  approximated  by  only  one  PDE  as 
shown  in  Eqs.  (35)  and  (46).  This  approximation,  however,  only 
results  in  a  little  loss  of  information  while  the  computational  time 
is  reduced  since  the  number  of  equations  that  are  solved  simulta¬ 
neously  is  reduced  by  a  factor  of  three. 

4.  Conclusion 

In  this  paper,  we  have  presented  a  method  for  simplifying 
rigorous  Li-ion  battery  models.  The  isothermal  model  used  for  sim¬ 
plification  is  based  on  the  first  principles  of  the  porous  electrode 
and  concentrated  solution  theories.  Besides  utilizing  the  nature  of 
the  battery  equations,  a  combination  of  several  techniques  such  as 
volume-averaging,  Galerkin’s  method,  and  curve-fitting  have  been 
used  to  achieve  an  approximate  model  that  can  simulate  in  mil¬ 
liseconds  without  a  significant  loss  in  accuracy  compared  to  the 
full-order  rigorous  model. 

The  paper  only  discusses  simulation  results  for  discharging 
behaviors.  However,  it  is  clear  that  the  method  presented  in  this 
paper  is  also  directly  applicable  for  charging  as  it  only  changes  the 
sign  of  the  applied  current  density. 
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